#include "AnaDef.h"
#include "AnaAliDef.h"
#include "dEdxGausN.h"

Int_t fgrID=-999;
Bool_t fklogx = 1;

//output in format of getGrSys
void getMCtrueFromPT()
{
  TH2D * htype = (TH2D*) gDirectory->Get("mcpdgVsPt");
  if(!htype){
    printf("no htype!! \n"); gDirectory->ls(); exit(1);
  }

  TH2D * hpars[]={0x0,0x0,0x0,0x0};
  for(Int_t ii=0; ii<4; ii++){
    hpars[ii]=(TH2D*) gDirectory->Get(Form("mcpar%dPt",ii));
    if(!hpars[ii]){
      printf("no mcpar %d !!\n", ii);  gDirectory->ls(); exit(1);
    }
  }

  TGraphAsymmErrors * grSys[4]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors * grStat[4]={0x0,0x0,0x0,0x0};
  TList *ll=new TList;
  for(Int_t ii=0; ii<4; ii++){
    grSys[ii]= new TGraphAsymmErrors;
    grSys[ii]->SetName(Form("grsys%d",ii));
    ll->Add(grSys[ii]);
    grStat[ii]= new TGraphAsymmErrors;
    grStat[ii]->SetName(Form("grstat%d",ii));
    ll->Add(grStat[ii]);
  }

  printf("check grDI %d\n", fgrID);
 
  //TF1 * fg=new TF1("fg","[0]*TMath::Gaus(x, [1], [2], 1)",0,500);

  for(Int_t ix=1; ix<=htype->GetNbinsX(); ix++){
    const Double_t xx = htype->GetXaxis()->GetBinCenter(ix);

    const Double_t nor = htype->Integral(ix,ix,1,4);

    //integral for MC in FF is not working due to cross-section weighting
    /*
//not for new IDFF
    if(!dEdxGausN::IsFF()){
      if(dEdxGausN::FilterIX(ix))
        continue;

      //nor only include 4 pars
      if(nor<dEdxGausN::PbinThres()){
        continue;
      }
    }
    */

    if(nor<EPSILON)
      continue;

    const Double_t ptmean = fklogx? (TMath::Power(10,xx)) : xx;
    //const Double_t cpion =  htype->GetBinContent(ix, dEdxGausN::kPION+1);

    for(Int_t iy = 1; iy <=4; iy++){
      const Int_t itype = iy-1;
      const Double_t cc = htype->GetBinContent(ix, iy);
      const Double_t ec = htype->GetBinError(ix,iy);

      //calculate others
      Double_t bb = 0;
      Double_t eb = 0;
      for(Int_t ib=1; ib<=4; ib++){
        if(ib==iy)
          continue;

        bb+= htype->GetBinContent(ix,ib);
        eb+= TMath::Power(htype->GetBinError(ix,ib),2);
      }
      eb = TMath::Sqrt(eb);

      /*
      //check consistency-->
      TH1D *h1d =hpars[itype]->ProjectionY(Form("pyx%dt%d",ix,itype),ix, ix);

      const Double_t check0 = h1d->Integral(0,100000);

      fg->SetParameters(TMath::Abs(h1d->GetEntries()*h1d->GetBinWidth(1)-10), h1d->GetMean(), h1d->GetRMS());
      Int_t kfit = h1d->Fit(fg,"LQ");
      const Double_t check1 = fg->GetParameter(0)/h1d->GetBinWidth(1);

      if(fabs(cc-check0)>1){
        printf("cc!=check0 %s %d %d:  %e -- %e %e\n", hpars[itype]->GetName(), ix, itype, TMath::Power(10,xx), cc, check0);exit(1);
      }
      if(fabs(cc-check1)>2 && cc>5){
        fg->SetParameters(TMath::Abs(h1d->GetEntries()*h1d->GetBinWidth(1)-10), h1d->GetMean(), h1d->GetRMS());
        h1d->Fit(fg,"L");

        TFile * ferr = new TFile(Form("ferr_%d_%d.root", ix, itype), "recreate");
        h1d->Write();
        fg->Write();
        ferr->Save();
        ferr->Close();
        printf("cc!=check1 %s %d %d:  %e -- %e %e -- %d %e %e %e\n", hpars[itype]->GetName(), ix, itype, TMath::Power(10,xx), cc, check1, kfit, fg->GetParameter(0), h1d->GetBinWidth(1), fg->GetParameter(0)/h1d->GetBinWidth(1));
        if(kfit==0)
          exit(1);
      }
      delete h1d;

      //<---
      */

      Double_t yy = -999;
      Double_t ey = -999;
      if(fgrID == 70 || fgrID==80 || fgrID==90){
        if(dEdxGausN::IsFF()){
          yy = cc/(cc+bb);

          ey = TMath::Sqrt( TMath::Power(cc*eb,2) + TMath::Power(bb*ec,2) )/TMath::Power(cc+bb,2);
          //ey = ec/(cc+bb); 
          printf("test itype %d pt %e cc %e ec %e bb %e eb %e nor %e -- yy %e ey %e\n", itype, ptmean, cc, ec, bb, eb, nor, yy, ey);
        }
        else{
          yy = cc/nor;
          ey = sqrt(cc)/nor;
        }
      }
      else if(fgrID == 23){
        printf("not implemented!!\n"); exit(1);
        /*
          const Double_t eta0 = dEdxGausN::Eta0();
          const Double_t eta1 = dEdxGausN::Eta1();
          
          printf("check eta %e %e\n", eta0, eta1);
        const Double_t etacorr = dEdxGausN::RatioEtaCorrection(itype,ptmean, eta0, eta1);

        printf("longtest etacorr %d %e -- %e\n", itype, ptmean, etacorr);

        yy = cc/cpion * etacorr;
        ey = yy*sqrt(1/cc+1/cpion) * etacorr;
        */
      }

      //printf("longtest %d %d -- %e -- %e %e %e \n", ix, iy, ptmean, cc, nor, yy);

      TGraphAsymmErrors *gr = grStat[itype];
      Int_t np = gr->GetN();
      gr->SetPoint(np, ptmean, yy);
      gr->SetPointError(np, 0,0,ey,ey);

      gr = grSys[itype];
      np = gr->GetN();
      gr->SetPoint(np, ptmean, yy);
    }
  }

  TFile *fout=new TFile(Form("fromPTmcfrac_%s_id%d.root", dEdxGausN::Opt(), fgrID),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();

  printf("done!\n");
}

int main(int argc, char * argv[])
{
  for(int ii =0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=4){
    printf("argc!=4 fin grID\n");return 1;
  }

  dEdxGausN::Ini();

  TFile::Open(argv[1]);
  fgrID = atoi(argv[2]);
  fklogx = atoi(argv[3]);

  if(fgrID!=70 && fgrID!=80 && fgrID!=90){
    printf("bad fgrID!! %d\n", fgrID);
    return 1;
  }

  getMCtrueFromPT();
}
